library(FactoMineR)
library(factoextra)
library("corrplot")
library(ggsci)
library(plotly)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(ComplexUpset)
library(pander)
library(knitr)
library(discoveR)
library(ade4)
library(adegraphics)
library(vegan)
library(MASS)
library(ComplexUpset)
library(RColorBrewer)

Dans tout le jeu de données, on retrouve :
- 9 stress biotiques (Biotrophic bacteria, Fungi, Nématodes, Oomycète, Rhodococcus, Stifenia, Necrotrophic bacteria, Virus et Other biotic).
- 9 stress abiotiques (Heavy metal, UV, Drought, Gamma, Nitrogen, Oxydative stress, Salt, Temperature et Other abiotic).

Seuls les stress biotiques sont étudiés ici.
Les fichiers utilisés ont deux colonnes informatives avec le stress appliqué sur l’échantillon en première colonne et le SWAP_ID pour pouvoir remonter aux informations de l’expérience précise.

Le fichier de chaque Gene Set est chargé à partir de la liste fournie dans le script. Les échantillons pour ce set sont séparés enuite pour les stress biotiques et abiotiques.

# Chemin vers dossier fichiers .txt et chargement
setwd("../docs_txt/")
data=read.table('Gene_Swap_NO_NA.dat',header=TRUE,sep='\t')

# Noms des colonnes sans les variables qualitatives
Gene=noquote(colnames(data[,c(-1,-2,-3)]))

#Noms des fichiers pour lecture dans boucle
GOSLIM=c("circadian_rythm.txt","abiotic.txt","biotic.txt","endogenous.txt",
         "external.txt","flower.txt","growth.txt","light.txt","photo.txt",
         "stress.txt")

matrix=cbind(rep(0,17341),Gene) #initialisation matrice

#Remplacer 0 par TRUE ou FAlSE pour chaque gene set
for(i in 1:10){
  datago=read.table(GOSLIM[i],header=TRUE,sep='\t')
  colgo=noquote(colnames(datago[,c(-1,-2,-3)]))
  test=match(Gene,colgo,nomatch =0)
  test[test!=0]<-"TRUE"
  test[test==0]<-"FALSE"
  matrix=cbind(matrix,test)
}

matrix=as.data.frame(matrix)
matrix=matrix[,-1]#enlever la première colonne de zéro d'initialisation
names(matrix)=c("Gene","Circadian","Abiotic","Biotic","Endogenous_stimulus",
       "External_Stimulus","Flower","Growth","Light","Photosynthesis",
       "Stress")
goslim=colnames(matrix)[2:11]

#Construction upset
(upset(matrix,goslim, name="GO SLIM", min_size=25,min_degree=1)
+ ggtitle("Tous les gènes (173441)"))

setwd("../docs_txt/")

#Variables importantes pour code, noms des fichiers pour lecture, nom des gene set pour entêtes personnalisées, différents stress..
files=c("circadian_rythm.txt","abiotic.txt","biotic.txt","endogenous.txt",
         "external.txt","flower.txt","growth.txt","light.txt","photo.txt",
         "stress.txt","random.txt")

name=c("Circadian rythm","Abiotic","Biotic","Endogenous stimulus",
       "External stimulus","Flower","Growth","Light","Photosynthesis",
       "Stress","Random")

biotic=c("BIOTROPHIC.BACTERIA","FUNGI","NECROTROPHIC.BACTERIA","NEMATODES","OOMYCETE","OTHER-BIOTIC","RHODOCOCCUS","STIFENIA","VIRUS")

#Boucle principale lecture Gene set
for(i in 1:11){
#Chargement Gene set 
  data=read.table(files[i],header=TRUE,sep='\t')
  #Conservation des échantillons seulement abiotiques
  data_biotic=data[data[,2] %in% biotic, ]
  row.names(data_biotic) <- data_biotic$Swap_ID
  
  #ACP avec données abiotiques
  res1=PCA(data_biotic[,c(-1,-3)],
           scale.unit=FALSE,
           graph=FALSE,
           quali.sup=1)
  fic.pca1<-dudi.pca(data_biotic[,c(-1,-2,-3)],
                     center=FALSE,
                     scale=FALSE,
                     scannf=FALSE)
  
  #BCA données abiotiques
  fic.bca1<-bca(fic.pca1,fac=as.factor(data_biotic$vec),scannf=FALSE)
  fic.tst<-randtest(fic.bca1)

#--------------------------------------------------------------------------------------------------#
                                        #Ecriture dans ouput HTML
#--------------------------------------------------------------------------------------------------#
  
  cat("\n")
  pander::pandoc.header(paste0("**Gene set : ",name[i]," (",dim(data)[2]-2," gènes)**\n"), level = 1)

#--------------------------------------------------------------------------------------------------#
                                                   #ACP
#--------------------------------------------------------------------------------------------------#  
  
  cat("\n")
  pander::pandoc.header(paste0("**ACP (SWAP_ID)**\n"), level = 2)
      
    g1<-fviz_pca_ind(res1,
             repel=TRUE,
             habillage = 1, # color by groups
             addEllipses = TRUE, # Concentration ellipses
             palette="Set1")
    print(g1)
      
    cat("\n")#Screeplot intertie
    pander::pandoc.header(paste0("*Inertie*\n"), level = 3)
    print(fviz_screeplot(res1, addlabels = TRUE, ylim = c(0, 50)))
    
    cat("\n")#Screeplot contribution individus et variables
    pander::pandoc.header(paste0("*Contribution*\n"), level =3)
    
    cat("\n")
    pander::pandoc.header(paste0("**Individus**\n"), level =4)
    print((fviz_contrib(res1, choice = "ind", axes = 1, top = 10)
        +
    fviz_contrib(res1, choice = "ind", axes = 2, top = 10)))
    
    cat("\n")
    pander::pandoc.header(paste0("**Variables**\n"), level =4)
    print((fviz_contrib(res1, choice = "var", axes = 1, top = 10)
        +
    fviz_contrib(res1, choice = "var", axes = 2, top = 10)))
    
    cat("\n")#Boxplot des meilleurs contributeurs
    pander::pandoc.header(paste0("**Distribution meilleurs contributeurs**\n"), level =4)
    #Récupération des labels des contributeurs via les plots factoextra
    #Découpage en 2 fois pour les deux dimensions    
    p=fviz_contrib(res1, choice = "var", axes = 1, top = 10)
    p2=fviz_contrib(res1, choice = "var", axes = 2, top = 10)
    contrib=p$data[order(p$data$contrib,decreasing=TRUE), ]
    contrib2=p2$data[order(p2$data$contrib,decreasing=TRUE), ]
    contrib=contrib[1:5,]
    contrib2=contrib2[1:5,]
    
    #Match entre noms des colonnes data_abiotic et les noms des meilleurs contributeurs    
    col_names=colnames(data_biotic)
    index=match(contrib$name,col_names)
    index=index-2
    index2=match(contrib2$name,col_names)
    index2=index2-2
  
      #Récupération des données de puces (sans colonnes qualitatives) pour les colonnes d'intérêt
    contrib_data=data_biotic[,c(-1,-2)]
    contrib_data=contrib_data[,index]
    contrib_data2=data_biotic[,c(-1,-2)]
    contrib_data2=contrib_data2[,index2]

    data_long=contrib_data%>%pivot_longer(
                         cols = c(1:5),
                         names_to = "Gene",
                         values_to = "Expression")
    
    g1<-ggplot(data_long)+
    geom_boxplot(aes(x=Gene, y = Expression))+
      theme(axis.text.x = element_text(angle = 60, hjust = 1))+
      ggtitle("Contributions Dim 1")
     
    
    
    data_long2=contrib_data2%>%pivot_longer(
                         cols = c(1:5),
                         names_to = "Gene",
                         values_to = "Expression")
    
    g2<-ggplot(data_long2)+
    geom_boxplot(aes(x=Gene, y = Expression))+
      theme(axis.text.x = element_text(angle = 60, hjust = 1))+
      ggtitle("Contributions Dim 2")
  
        #Détermination des limites en fonction des valeurs des deux dimensions pour rendre comparaison plus facile
    g1<-g1+ylim(min(data_long$Expression,data_long2$Expression),
                max(data_long$Expression,data_long2$Expression))
    g2<-g2+ylim(min(data_long$Expression,data_long2$Expression),
                max(data_long$Expression,data_long2$Expression))
         print(g1+g2)

#--------------------------------------------------------------------------------------------------#
                                                #BCA
#--------------------------------------------------------------------------------------------------#  
         
      cat("\n")
      pander::pandoc.header(paste0("**BCA**\n"), level = 2)
      cat("\n")
      pander::pandoc.header(paste0("*Toutes positions*\n"), level = 3)
      cat("\n")
      
      .z<-paste("BCA sur position \n ratio= ",round(fic.bca1$ratio,2)," pval= ",signif(fic.tst$pvalue,digits=2))

      s.class(fic.bca1$ls,
              fac=as.factor(data_biotic$vec),
              col=rainbow(9),
              psub=list(text=.z,cex=1,position="topleft"),
              plabels.cex=0,
              key.cex=0.8,
              key.text.cex=0.8)

      cat("\n")
      pander::pandoc.header(paste0("*Par position*\n"), level = 3)

      #Plot pour chaque stress     
      plotlist=list()
      for(j in 1:9){
        s.class(fic.bca1$ls,
                fac=as.factor(data_biotic$vec),
                psub=list(text=paste("BCA \n", biotic[j]),cex=1,position="topleft"),
                plabels.cex=0,
                key.cex=0.4,
                key.text.cex=0.4,
                col=c(rep(grey(0.95),j-1),brewer.pal(n=9,name="Set1")[j],rep(grey(0.95),9-j)),
                plot=FALSE,
                chullSize=0,
                ellipseSize=0)->plotlist[j]
      }

      ADEgS(plotlist,layout=c(3,3))
}

Gene set : Circadian rythm (58 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Abiotic (662 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Biotic (421 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Endogenous stimulus (524 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : External stimulus (553 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Flower (187 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Growth (223 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Light (293 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Photosynthesis (76 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Stress (1178 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position

Gene set : Random (51 gènes)

ACP (SWAP_ID)

Inertie

Contribution

Individus

Variables

Distribution meilleurs contributeurs

BCA

Toutes positions

Par position